QCD saturation at the LHC: comparisons of models to p+p and A+A data and 

predictions for p+Pb collisions 
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In a previous paper (larXiv:1011.1895"|l . we showed that saturation models, constrained by e+p 
HERA data on inclusive and diffractive cross-sections, are in good agreement with p+p data at 
LHC in the soft sector. Particularly impressive was the agreement of saturation models with the 
multiplicity distribution as a function of n^h.- In this paper, we extend these studies further and 
consider the agreement of these models with data on bulk distributions in A+A collisions. We 
compare our results to data on central and forward particle production in d+Au collisions at RHIC 
and make predictions for inclusive distributions in p+Pb collisions at the LHC. 



I. INTRODUCTION 

In a previous paper [l| , we discussed the computation of inclusive distributions in p+p collisions at RHIC and the 
LHC within the framework of k± factorization, wherein the unintegrated gluon distributions were determined from 
fits to small X HERA data on inclusive, diffractive and exclusive final states. The key ingredient in these fits is the 
dipole cross-section, which, to leading logarithmic accuracy, can be defined as 
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where U(h± ± is a Wilson line in the fundamental representation representing the interaction between a quark 
and the color fields of the target. In the Color Glass Condensate (CGC) framework 0,0] of gluon saturation [J, the 
average {■ ■ ■)x is an average over these color fields; the energy dependence of the correlator as a function of x (or 
the rapidity Y = ln(l/a;)) is given by the JIMWLK equation 15|. In the large Nc limit, the equation for the energy 
evolution of this correlator is the Balitsky-Kovchegov (BK) equation Q- 

We note however that neither JIMWLK nor BK is at present equipped to deal well with the impact parameter 
dependence of the dipole cross-section; the dipole cross-section in this formalism is taken in eq. ^ to be independent 
of the impact parameter. Another limitation of this framework is that the full next-to-leading logarithmic (NLL) 
expressions are not yet available; at the NLL level, only running coupling corrections to the leading log kernel have 
been considered in phenomenological applications. With these limitations in mind, we considered in ref. saturation 
models of the dipole cross-section with the common criteria that their parameters be strongly constrained by fits to 
the HERA data. The saturation models considered^ included the IP-Sat model the b-CGC model [sl-fioj and the 
rcBK model [llj. All of these models provide good fits to the HERA data^. 

In hadron-hadron collisions, one can derive at leading order the expression (l5j 
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This equation is a generalization of the well known k±_ factorization expression for inclusive gluon production 
to include the impact parameter dependence of the unintegrated gluon distributions. Here Cp = {N^ — i)/'2Nc is 
the Casimir for the fundamental representation. Using a relation between quark and gluon dipole amplitudes strictly 
valid in the large Nc limit, the unintegrated gluon distribution in either of the two protons can be expressed in terms 
of the corresponding dipole cross-section measured in DIS as [l7j 
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^ For more details, sec the discussion in section II of ref. 

^ The rcBK model has only been compared to inclusive HERA data [T3l . We note however that this comparison is to the combined 
Hl-ZEUS inclusive data [l3{ . in contrast to the IP-Sat and b-CGC models, which were fit only to the older ZEUS data [ij- 
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Thus the impact parameter dependent dipole cross-section determined from HERA data can be used to compute 
single inclusive gluon distributions in proton-proton collisions. Since this is a leading order computation, the overall 
normalization is not constrained and is determined from data as we shall describe below. For the integrated mul- 
tiplicities, there is a logarithmic infrared divergence that can be regulated by introducing a mass term as discussed 
in our previous paper. We should mention that solutions of Yang-Mills equations that treat the infrared behavior 
properly give infrared finite distributions [l8l - [20l | . For a nice comparison of theoretical errors in various k± factorized 
approximations to the full classical Yang-Mills results, see ref. [21[ . 

In rcf. P], we used eq. ^ combined with eqs. ^ and ([T]) to compute the rapidity and p± distributions for the 
models discussed for central rapidities in p-l-p collisions at RHIC energies all the way to the highest available LHC 
energies. A typical feature of the rapidity and p± distributions in these models was that the agreement with data 
improved with increasing energy; this should be the case because saturation effects are increasingly important at 
higher energies. We also observed that the saturation models, in particular the IP-Sat model, gave excellent fits to 
the multiplicity distribution P{nch.) as a function of rich, for a wide range of energies. In this work, we will first 
briefly consider p+p collisions again before discussing A-l-A collisions and p/d+A collisions. In the latter case, we 
will present predictions for a future p-(-Pb run at the LHC. 

II. RESULTS FOR P+P COLLISIONS 

The probability distribution for producing n particles is 

P{n) = I d2b^^pNB(-(b^),fc(b^)), (4) 

where, in the Glasma fiux tube framework [2^ . we can derive the negative binomial distribution (23l | 

pNB/- _ T{k + n) n^k^ 
" ^ ' ^ " T{k)T{n + l){n + ky-+^ ' ^ ' 

with the parameter k defined specifically to be 

Mb.)^c ''"'-;5"\ (0) 

with C=0.155 obtained from a fit to p+p multiplicity distribution. Here 



where, motivated by CGC computations on dilute-dense collisions [2J], we choose Qs in the overlap area of the two 
hadrons to be Qs(s±, bj^) — viAn.{Qs{s±),Qs{si. — bj,)}. Also, n is the average multiplicity at a given impact 
parameter in the saturation model. Finally, as previously, we use 
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In fig. ([T]), we show results for the p+p multiplicity distribution plotted for ^/s up to 7 TeV for \r]\ <0.5. This 
extends the comparison in ref. [l| to the CMS data at 7 TeV, with good agreement to rich. = 60, nearly twice the 
range considered previously. The parameter ^=0.155 is extracted from a fit to data at 0.9 TeV and used for all other 
energies^. It is the only free parameter in our fit; however, this quantity was recently computed non-perturbatively ab 
initio by solving Yang-Mills equations numerically [26l | for two gluon correlations from gauge fields generated in the 
collision of two dense color sources. The results of the numerical computation vary depending on parameter choices in 
the range ( ^ 0.3-1.5-the lower end of this range is therefore a factor of two larger than the best fit value. Given the 
many uncertainties in the computation, the agreement is quite good, keeping in mind that nothing a priori prevents 



^ To avoid confusion, the value for (J quoted in ref. [l[ is for in the fundamental representation; fc(bx) is identical to that discussed 
here. Our estimations in the adjoint representation are shown in ref. |25|| . 
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C from being orders of magnitude different. The good agreement of this framework with the LHC data over several 
decades, taken at face value, leads us to conjecture'* that fluctuations in the number of produced gluons for a fixed 
distribution of hot spots gives a larger contribution to the multiplicity distribution than fluctuations in the distribution 
of hot spots themselves. The latter is of course only treated at the mean field level here in contrast to "pomeron loop" 
contributions [l^l -suggesting perhaps that the latter are suppressed. Because Glasma flux tubes generate long range 
rapidity correlations, and can explain the distribution of high multiplicity events, our result provides further evidence 
corroborating computations in this framework [28l . [29j that suggest Glasma flux tubes generate the near side ridge 
seen in high multiplicity events by the CMS collaboration [s^l ■ 
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FIG. 1: Multiplicity distribution for p+p collisions in the IP-Sat model compared to data from UA5, CMS and ALICE [31l43i 



We now turn to transverse momentum distributions for charged hadrons and tt s, that in p+p collisions are 
computed from the expression 
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where^ Dg-^h{z, p^) is chosen to be 6.05 z-O-'^^''^^ - zf-^^. In ref. |H, we had considered the transverse momentum 
distributions only at mid-rapidity, where agreement of the IP-Sat and rcBK models was not particularly good. This 
agreement improved significantly at the higher LHC energies. The explanation provided there was that this better 
agreement is a consequence of smaller x values being probed at the higher energies. In fig. we show the results of 
the rc BK and IP-Sat models for forward rapidities for RHIC p+p collision at 200 GeV. 

For both IP-Sat and rcBK we have used an overall normalization^ extracted from energy dependence of the single 
inclusive multiplicity of the form: A/{'Kb'^^y^ with &max = bQ + C ln{y/s). This form absorbs the uncertainties in the 
inelastic cross-section and higher order effects (iiT- factors). For the IP-Sat (rcBK) model, one finds A = 0.23(6.15), 
bo = 5.77(5.14) and C = 0.32(0.76) using the mass term m = 0.4 GeV by fitting data points at 77 = over the 
range of energy shown in fig. [31 For the rc-BK model the constant term A absorbs the prefactors ^ of Eq. [2] which 
includes unknown overlap area and other terms that cannot be separated from the " K factor " . The results are 
rather insensitive to the infrared cut-off to. However, one finds a ^ 10% variation of the normalization when the 



* To confirm this conjecture, one would need to demonstrate that the results are valid for instance in the rcBK framework. However, 

because the impact parameter dependence of inclusive distributions in this model is rather simplistic and not constrained by data, this 

is difficult to confirm meaningfully at present. 
^ The functional form quoted is for 7r+ -|- tt" ; we assume that the likelihood to fragment is identical for either charge and is equal to that 

for TT^. For other charged hadrons, the functional form is assumed to be identical with the normalization constrained by the momentum 

sum rule. 

^ The overall normalization for minimum bias p+p is slightly different from our previous paper Q due to one less parameter (Aq) in the 
large x extrapolation which is set to zero here. 

^ Eq.[2]for min-bias p+p collision the rc-BK model normalization constant includes the constant prcfactor ( ^^'^f -, — „ ' i^',^„'> , ) ; here 

Ra, Rb corresponds to the radii of two protons and Sa,b is the overlap area. 
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constants are extracted by a) fitting tlic full pscudo rapidity at RHIC energy, b) considering data points from different 
experiments. This variation, along with the numerical uncertainties, contributes to the gray bands shown in fig.[21 We 
see that the agreement at forward rapidities is significantly better than our previous comparison to the mid-rapidity 
distribution; this result provides a good benchmark for computing RpA at RHIC and in predictions of the same for 
p+Pb collisions at the LHC. 
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FIG. 2: Transverse momentum distributions at forward rapidities in rcBK and IP-Sat models compared to STAR [sj] and 
BRAHMS [3^ data. The gray bands show the uncertainty in the determination of the normalization constant. 



III. RESULTS FOR A+A COLLISIONS 



For a large nucleus, in the IP-Sat model, we can approximate the dipole-nucleus cross section to be 
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where crdip(r_L, a;)^ is obtained from integrating the dipole-proton cross section in cq. ([T]) over the impact parameter 
distribution in the proton. This form of the dipole-nucleus cross-section was shown previously to give reasonable 
fits to the limited available fixed target e-f A inclusive data. The initial conditions for rcBK evolution for a nucleus 
were similarly fixed by comparisons to the e-|-A data [s^] . 

Substituting the expression for the dipole-nucleus cross-section in eqs. ([2]), and likewise the latter in one can 
compute the nuclear multiplicity distributions. The infrared divergence in the multiplicity distribution is regulated in 
exactly the same was as was the case for the p-l-p multiplicity distribution, by replacing p± by mj^ = ■\/p5_ + fn'^, with 
m varied between 0.2-0.4 GeV. Wherever we have considered fixed coupling, we have used a5=0.2; for the running 
coupling case, we run as with the scale Qs ~ max. {Qsi^i, s±), Qs(x2, s± ~ b^)}. 

Fig. ([3]) shows the energy dependence of average multiplicity for most central Au+Au collision for fixed and running 
coupling in the IP-Sat model. The number of participants^ at a given impact parameter is determined from the Glauber 
relation 13811 



iVpart(b±) ^ A j Ta{s^) {1 - [1 - TBis^-hi_)aNnf}d^s^ + B j Tb{si_) {l - [1 - TaIs^- hi_)<j^^]^} d^s 



The results shown in fig. ^ are for 0-6% centrality, which corresponds to a median ~ 12.2 GeV^^; we compute 
'^^°'d!)°^^ ^'^d ^part(bi) for this median value. We observe that a fairly good agreement with data is obtained for the 
infrared cut-off given by m = 0.4 GcV. The prescription for the running coupling gives a variation that corresponds 
to a 20% uncertainty at lower energies, which decreases significantly at higher energies. 



In this expression, ctnn ~ 62 mb for 2.76 TeV and (Tnn ^ 41 mb at 200 GeV. 
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FIG. 3: Left: Energy dependance of the multiplicity per participant in the IP-Sat model for p+p and A+A collisions. For the 
A+A case, the calculation is done for the — 6% centrality. Right: same plot for A+A with fixed (solid) and running (dashed) 
coupling. Data points for p+p are from ref. [39l - [43 | and for A+A from ref. [13, \EM 



We next consider the centrality dependence of the multiplicity at RHIC {y^s = 200 GeV) and LHC {y^ =2.76 
TeV) in the IP-Sat model. While the agreement of the model with data shown in fig. (U) is reasonably good for the 
most central collisions, a systematic deviation is seen for lower centralities, and the model underpredicts the data. 
While within the range of the theoretical uncertainties outlined thus far, this systematic discrepancy leaves significant 
room for final state entropy production, which is expected to be more significant for more peripheral collisions. See 
for instance refs. [i^ that estimate the amount of entropy production. As the right plot of fig. (jlj shows, running 
coupling effects are less important for the most central collisions but introduce significant uncertainties relative to the 
fixed coupling results for more peripheral collisions. 
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FIG. 4: Centrality dependance of the inclusive multiplicity in the IP-Sat model compared to RHIC [47j and LHC [i^ data. 
Left: (fixed coupling) 200 GeV values for both data and model are multiplied by a factor 2.08. Right: Same plot comparing 
running (solid curve) and fixed coupling (dashed curve) results in the IP-Sat model. 



Fig. ([5]) shows the pseudo-rapidity distributions in the IP-Sat model compared to data for Au+Au collisions at 
200 GeV (PHOBOS) and Pb+Pb collisions at 2.76 TeV (ALICE and CMS). Firstly, one sees that the results are 
sensitive to the infrared cut-off, with improved agreement seen for m = 0.4 GeV. Further, the rapidity distributions 
are sensitive to the extrapolation of the model to larger x > 0.01 values. We also note that a significantly better fit 
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to the data at higher energies is obtained by including running couphng effects. We now employ cq. (|4]) to compute 
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FIG. 5: Pseudo-rapidity distribution from kx-factorization formula. Left: fixed coupling results for 200 GeV and 2.76 TeV. 
Right: Same plot(m=0.4) at 2.76 TeV with running (solid curve) and fixed coupling(dashed curve) in kx-factorization formula. 
Data points are from ref. [49l . [sot 



the multiplicity distribution in A+A collisions. While eq. ^ is computed identically to the p+p case, we need to 
determine the impact parameter distribution differently from the prescription used for the p+p case in cq. ([T])- The 
expression 



1 - (1 - cnnTab) 



AB 



(10) 



gives a better description of the impact parameter distribution in A-|-A collisions^. As in the computation of iVpart, 
ctnn ^ 62 (41) mb for 2.76 (0.2) TeV, is standard and not varied. The saturation scale in this computation is 
determined at the median value of the impact parameter b_L™'"^' = 15 GeV~^. With these assumptions, the only 
parameters in computing P{n) are m and C, the parameter controlling the width of the multiplicity distributions. In 
fig. (HI) (left) we see that the multiplicity distributions are insensitive to variations in m. Fig. 1^ (right) shows the 
result of varying ( = 0.01-1. Interestingly, wc find that the best fit is found for the value of C = 0.155 that also gives 
the best fit to the p + p data. 

For the A-l-A collision data discussed here, we presented results for a) the energy dependence of central (0-6%) 
A-l-A collisions, b) the centrality dependence at two different energies, c) the rapidity dependence (at RHIC for 0-6% 
centrality and at LHC for 0-5% centrality), d) the multiplicity distributions (which are minimum bias). Unlike p+p 
collisions, for A+A collisions, we do not use an energy dependent normalization (because the nuclear size does not 
grow appreciably with energy unlike the proton) but a constant factor that is fit to the single inclusive distribution 
at one energy at = and used to fit these four data sets a) -d). For the IP-Sat model (the only model considered 
here), the value of this normalization (K-iactor) is 0.07 for the mass term m=0.4 



^ For A+A collisions, the published data |43 | points are uncorrected requiring an additional parameter in contrast to the p+p case. The 
average multiplicity is 



n(bx) = C„ 



dN{bj_) 



drj 



(11) 



|i7l<0-5 



where the pre-factor Cm is the additional parameter specific to the multiplicity distribution and is tuned to provide a good fit to the 
uncorrected multiplicity distribution. 
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FIG. 6: Left:Multiplicity distribution for Au+Au collisions in the IP-Sat model compared to uncorrected data(liistogram [isl ]) 
for different values of m. Both data and model plots are normalized for better comparison. Right:Multiplicity distribution for 
Au+Au collisions at 200 GeV and its sensitivity to the non perturbative constant 



IV. RESULTS FOR p+A COLLISIONS 



We computed the min-bias average multiplicity at mid-rapidity for the IP-Sat and rc BK models for p-|-A collisions 
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The data is normalized to the PHOBOS 200 GeV d-|-Au data [51[. The energy dependence of the average multiplicity 
is shown in fig. ([7]), the band corresponds to the variation of m in the range of 0.2-0.4 GeV. Both the models give a 
comparable energy dependence, with the IP-Sat model giving a slightly higher multiplicity at the highest energies. The 




FIG. 7: Energy dependence of the minimum bias single inclusive multiplicity at ry = in p+A collisions from fcx-factorized 
unintegrated distributions determined in the rcBK and IP-Sat models. The distribution is normalized with respect to the 
PHOBOS d-|-Au data [sj] for 200 GeV. The band represents the uncertainty in the calculation due to the variation of the 
mass term in the range of 0.2-0.4 GeV. 



rapidity distributions for RHIC energies in the two models and predictions for LHC energies are shown in fig. ([S]). The 
models agree with the RHIC data with an accuracy of sa 10%, which is within the theoretical systematic uncertainty, 



The non-perturbative scales at the initial rapidity for the rc-BK model [i[ used here are Q^q ^=0.09 GeV^ and ^=0.44 GeV'^. 
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which we shall discuss further shortly. 
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FIG. 8: Pseudo-rapidity distribution for minimum-bias p-l-A collision at RHIC and LHC energies. Prediction from rcBK shown 
for mass term m=0.2 and 0.4 GeV. Data points are from ref. [sil, [s^l 



Fig. ^ shows the transverse momentum distribution compared to BRAHMS [35| and STAR [34| 200 GeV data 
for and tt^ at forward rapidities. For the h~ case, we have included 15% isospin correction in the normalization 
constant, to be discussed further below. Predictions for the p + A transverse momentum distributions for charged 
hadrons at 77 = for LHC energies with this fixed normalization are shown in fig. 
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FIG. 9: Transverse momentum distribution at forward rapidity at the highest RHIC energy compared to STAR [3J] and 
BRAHMS [s^ data. The gray bands show the uncertainty in determination of normalization constant from various sources. 



We compute the nuclear modification factor anticipated in p + A collisions at the LHC in the IP-Sat and rcBK 
models. RpA in both models, for y/s = 4.4 TeV/nucleon, approaches unity at p_L ~ 5 GeV. For y/s = 8.8 TeV/nuclcon, 
the suppression persists to higher p±. The slope of RpA however appears quite different in the two models. 

Regarding the overall normalization, for d-f Au we have data available at only one energy. As noted, data and 
predictions include a) the energy dependence, b) the rapidity dependence at RHIC and LHC, and the the single 
inclusive pt distributions. 

For the IP-Sat model the form of the normalization is A/lirb'^^y^) where, A = 0.25 with mass term m = 0.4 GeV. 
Unlike p + p here 6max = 9.5 fm is a fixed number which is the maximum range of impact parameter to obtain 
minimum bias distribution and does not change with energy. For the rc-BK model a single normalization constant 
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FIG. 10: Transverse momentum distribution at mid-rapidity for minimum-bias p-l-p and p-l-A collisions.. 
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FIG. 11: Nuclear modification factor for charged hadrons at mid-rapidity for p-|-A collisions at LHC energies. 

A = 0.032 (with m = 0.4 GcV) absorbs all the constant pre factors of eql2l These normalization constants are 
obtained from a fit to the PHOBOS pseudo rapidity distribution; one obtains an ~ 8% higher value of A when the 
fit is performed only to the data points at 77 = 0. Using the BRAHMS data for normalization also gives higher 
values for A which, along with other numerical uncertainties, contributes to the bands shown in fig. [9] The significant 
difference in A for IP-Sat and rcBK for d-|-Au collisions is because the the area of overlap and other terms in the 
pre-factor of kt-factorization are absorbed in A for rc-BK and cannot be separated from the "K factor" . For the 
IP-Sat model, eq. ((3)) includes those factors and A is of the order of 1. This apparent difference in the two models 
doesn't affect any of our final results since same normalization is consistently used everywhere. Note also that in 
conversion from d-|-Au to p-|-A numbers we have used an additional factor of 1.6/2 in the normalization, which is 
standard in such conversions in the literature. In computing RpA, the result depends on A^coii, which is sensitive to the 
proton inelastic cross-section. Since 6max,proton grows with energy, one finds for the energies -y/s = 4.4, 8.8 TeV that 
ratio 6max.proton(8-8TeV)/6^-,ax,ppton(4-4TeV) 7Vcoii(8.8TeV)/iVcoii(4.4TeV), a result consistent with expectations of 
TVcoU from Glauber approaches [53|; numbers quoted are in agreement with our ratio to 5%. 



eqlUfor min-bias p+A collision the rc-BK model normalization constant includes the pre factor ( ,e -, — „,f;'^„.) , ): here Rr,, Ra 

corresponds to the radii of the proton and nucleus and Sp,A is the overlap area. 
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V. CAVEATS AND COMPARISONS 



We would be remiss not to discuss the known sources of theoretical uncertainties in our framework; a brief discussion 
was also presented in ref. [l|. One that we alluded to previously is the kj_ factorization framework. Corrections to 
these come both from additional multiple scattering contributions in the dense projectile-dense target limit and from 
higher order contributions [13, ; We estimate differences with the classical Yang-Mills multiple scattering effects to 
be no larger than other theoretical uncertainties and is accounted for an approximately 15% difference in normalization 
between the integrated multiplicity and the p± distributions. The higher order contributions are constrained by a 
normalization constant that's fit at a given energy. We have also shown in fig. ([S]) the effect of running coupling on 
multiplicity distributions. Another source of theoretical uncertainty arises from logarithmic sensitivity to the infrared 
cut-off seen in the single inclusive multiplicity (but not in the multiplicity distribution); the additional (classical 
Yang-Mills) scattering effect we alluded to also makes the distributions infrared finite. Even though the Yang-MIlls 
contribution is understood, it is at present cumbersome to include it in a "global" analysis. 

In the IP-Sat model, as discussed previously in ref. [l[ , there is a sensitivity to the different data sets for the dipole 
cross-section that gave best fits to the then extant HERA data. We noted in ref. [ij that this sensitivity decreased 
with increasing energy. The sensitivity to the different fits to the giuon distribution is particularly severe at large p± 
where small differences are accentuated because of the steeply dropping cross-section. At high p± and large rapidities, 
the unintegrated distributions are sensitive to a; > xo{= 0.01). Guided by quark counting rule prescriptions, we choose 
(f>{x > xo) = (f>{xo){l — a;)'^/(l — xq)^, as in eq. (20) of ref. [l|, except with the parameter A there set to zero. Further, 
as noted previously, the IP-Sat fits in ref. ^ were performed prior to the combined ZEUS-Hl data . One expects 
a revised fit to have an impact on the numbers extracted in ref. Finally, there is theoretical uncertainty arising 
from our imperfect knowledge of fragmentation functions. 

Given the many uncertainties, one has to think of these predictions at this stage as having considerable elasticity of 
20%-50% variability depending on the quantity studied. Only qualitative differences in measured distributions from 
theoretical projections (such as RpA > 1 at low p± at the LHC) can rule out models at present. Else, one has to check 
whether variations in parameters without adding new ones can accommodate the data. This is no different in spirit 
from global fits in perturbative QCD. Given that we have not tried to fine tune parameters, the agreement with the 
data presented here (in combination to the results in ref. [l|) is quite good considering that the dipole approach also 
gives a good description of small x HERA data. 

We will now briefly compare our results with some recent related approaches. In ref. [56l |. results for the pseudo- 
rapidity distributions in p+p and A+A collisions at RHIC and the LHC are presented (based on previous work in 
ref. |57|) and predictions made for p-f A collisions. The results shown are primarily for the b-CGC model which we 
also considered previously. While the b-CGC model does well for the pseudo-rapidity distributions, it does less well 
relative to the IP-Sat model for the multiplicity distribution P{n); as fig. ([T|) indicates, the latter works quite well. 
When it comes to nuclei, ref. [56j assumes that an additional gluon cascade mechanism to explain the faster energy 
dependence of the charged particle multiplicity. However, as pointed out in ref. (ssj . the difference in the energy 
dependence of the charged particle multiplicity in A-l-A collisions relative to p-|-p collisions is economically explained 
by the fact that the larger saturation scales in nuclei lead to more phase space for bremsstrahlung; while this may be 
conceptually similar to the cascade posited in ref. [56j . no further modification of the formalism is necessary in the 
former approach. Indeed, we see in fig. that a good agreement is obtained for both p-l-p and A-l-A data in the 
IP-Sat model. A similar quality of agreement"'^^ is seen in the rcBK model (combining results in refs. [l| and (soj). 
In ref. [g^], results in the KLN model for the pseudo-rapidity distributions in p+p, d+A and A+A collisions at RHIC 
are presented, as are results for the LHC in p+p and A+A collisions, and predictions are made for the same in p+A 
collisions at the LHC. Similarly, results are presented for the P{n) in p+p collisions and predictions of the same made 
for p+A collisions. The KLN model gives good results for the pseudo-rapidity distributions at RHIC energies and in 
p+p and A+A distributions at LHC energies, which are similar to the models considered here. While the multiplicity 
distribution P{n) is motivated similarly to our discussion, the expression in ref. [60j does not depend on impact 
parameter; it will be interesting to see if a similar quality of agreement is obtained for higher n as in our fig. ([T]). The 
KLN model gives a better agreement to the centrality dependence of the charged particle multiplicity when compared 
to fig. (HI ; our result however leaves open the possibility that there may be significant entropy generation in the final 
state in more peripheral collisions. We also note that the KLN model does not at present allow for a global fit because 



it is not constrained by HERA data on e+p collisions [61 



Multiplicity distributions are sensitive only to p_L < 1 GeV. For p± distributions, we first compare our results for 



Small differences in the predictions there for pA pseudo-rapidity distributions with those presented here can be attributed primarily to 
slight different initial values for the saturation scale. 
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RpA to those in ref. [64I. Only results for 8.8 TeV/nucleon are presented there, and for forward rapidities. Even 
so, we see that the suppression is considerably less in our framework. A significant fact is that in RpA, one has to 
account for a growth in the proton minimum bias cross-section with energy, which tends to suppress the denominator 
and lessen the suppression from what it would be otherwise. We found this effect to be quite significant and modify 



plots from those similar to rcf. [62| to what we present here. We expect the same effect to impact the result (in 



fig. (5) of ref. [63| using a "hybrid" formalism [6J] corrected by the inelastic contributions introduced in ref. [65[. The 
"hybrid" approach may be more suitable than our approach when distributions arc sensitive to large x > 0.01 and 
high ; significant recent progress [55] in that approach will allow us to compute systematically corrections beyond 
the k± factorization framework. Finally, we observe that RpA in the IP-Sat model approaches unity more rapidly for 
= 4.4 TeV/nucleon than for the leading twist shadowing evolved EPS08 curve [H, |63| plotted in ref. [HI], while 
the rcBK curve has a similar behavior to the EPS08 curve. This is not too surprising because shadowing in the IP-Sat 
model is from higher twist contributions that go away quickly at large Q^; in contrast, there is a significant leading 
twist window in the rcBK model which may explain the qualitative similarity of the two results in this kinematic 
window. 



Acknowledgements 

R.V is supported by the US Department of Energy under DOE Contract No.DE-AC02-98CH10886. We thank 
Subhasis Chattopadhyay and Adrian Dumitru for very valuable conversations. Wc thank T. Lappi for bringing a typo 
in the earlier version of this paper to our attention. 



P. Tribedy, R. Venugopalan, Nucl. Phys. A850, 136-156 (2011) 
L.D. McLerran, R. Venugopalan, Phys. Rev. D 49, 2233 (1994); ibid. 49, 3352 (1994); ibid. 50, 2225 (1994) 
E. lancu, R. Venugopalan, Quark Gluon Plasma 3, Eds. R.C. Hwa, X.N. Wang, World Scientific 
H. Weigert, Prog. Part. Nucl. Phys. 55, 461 (2005); F. Gelis, E. lancu, J. Jalilian-Marian, R. Venugopalan, 



hep-ph/0303204 



arXiv:1002.0333 



L.V. Gribov, E.M. Levin, M.G. Ryskin, Phys. Rept. 100, 1 (1983); A.H. Mueller, J-W. Qiu, Nucl. Phys. B 268, 427 
(1986). 

Jalilian-Marian, A. Kovner, A. Leonidov, H. Weigert, Nucl. Phys. B 504, 415 (1997); ibid., Phys. Rev. D 59, 014014 
(1999); E. lancu, A. Leonidov, L.D. McLerran, Nucl. Phys. A 692, 583 (2001); E. Ferreiro, E. lancu, A. Leonidov, L.D. 
McLerran, Nucl. Phys. A 703, 489 (2002). 

I. Balitsky, Nucl. Phys. B 463, 99 (1996); Yu.V. Kovchegov, Phys. Rev. D 61, 074018 (2000). 
H. Kowalski, D. Teaney, Phys. Rev. D 68, 114005 (2003). 

E. lancu, K. Itakura, S. Munier, Phys. Lett. B 590, 199 (2004). 

H. Kowalski, L. Motyka, G. Watt, Phys. Rev. D 74, 074016 (2006). 

G. Watt, H. Kowalski, Phys. Rev. D 78, 014016 (2008). 

J. L. Albacete, Y. V. Kovchegov, Phys. Rev. D 75, 125021 (2007). 

J. L. Albacete, N. Armesto, J. G. Milhano, P. Quiroga- Arias and C. A. Salgado, Eur. Phys. J. C 71, 1705 (2011). 

F. D. Aaron et al. [HI and ZEUS Collaboration], JHEP 1001, 109 (2010). 
S. Chekanov et al., [ZEUS Collaboration] Eur. Phys. J. C 21 (2001) 443. 
J. P. Blaizot, F. Gelis, R. Venugopalan, Nucl. Phys. A 743,57 (2004). 
M.A. Braun, Phys. Lett. B 483, 105 (2000). 

F. Gelis, A. M. Stasto, R. Venugopalan, Eur. Phys. J. C 48, 489 (2006). 
A. Krasnitz and R. Venugopalan, Phys. Rev. Lett. 84, 4309 (2000). 
A. Krasnitz, Y. Nara and R. Venugopalan, Nucl. Phys. A 727, 427 (2003). 
T. Lappi, Phys. Rev. C 67, 054903 (2003). 

J. -P. Blaizot, T. Lappi and Y. Mehtar-Tani, Nucl. Phys. A 846, 63 (2010). 

A. Dumitru, F. Gelis, L. McLerran, R. Venugopalan, Nucl. Phys. A 810, 91 (2008). 

F. Gelis, T. Lappi, L. McLerran, Nucl. Phys. A828 (2009) 149. 

A. Dumitru and L. D. McLerran , Nucl. Phys. A 700, 492 (2002). 

P. Tribedy and R. Venugopalan. larXiv:1101.5922 [hep-ph]. 

T. Lappi, S. Srednyak, R. Venugo palan, JHEP 1001 066 (2010). 

A. Kovner and M. Lublinsky, arXi v:1109.0347l [hep-ph]. 

A. Dumitru, K. Dusling, F. Gelis, J. Jahhan-Marian, T. Lappi and R. Venugopalan, Phys. Lett. B 697, 21 (2011). 

K. Dusling and R. Venugopalan, arXiv: 120 1.2658 [hep-ph]. 

V. Khachatryan et al. [CMS Collaboration], JHEP 1009, 091 (2010). 

UA5 Collaboration, Z. Phys. C43 (1989) 357. 

ALICE Collaboration, Eur.Phys.J.C68:345-354,2010. 



12 



V. Khachatryan et al. [CMS Collaboration], JHEP 1101, 079 (2011). 
J. Adams et al. [STAR Collaboration], Phys. Rev. Lett. 97, 152302 (2006). 
I. Arsene et al. [BRAHMS Collaboration], Phys. Rev. Lett. 93, 242303 (2004). 
H. Kowalski, T. Lappi, R. Venugopalan, Phys. Rev. Lett. 100, 022303 (2008). 
K. Dusling, F. Gehs, T. Lappi and R. Venugopalan, Nucl. Phys. A 836, 159 (2010). 
D. Kharzeev and M. Nardi, Phys. Lett. B 507, 121 (2001). 
K. Aamodt et al. [ALICE Collaboration], Eur. Phys. J. C 68, 89 (2010). 
V. Khachatryan et al. [CMS Collaboration], JHEP 1002, 041 (2010). 
K. Aamodt et al. [ALICE Collaboration], Eur. Phys. J. C 68, 345 (2010). 



UA5 collaboration, G. J. Alner et. al, \Z. Phys. C33 (1986) T| UAl collaboration C. Albajar et. al. 



Nucl. Phys. B335 (1990) 261, CDF collaboration, F. Abe et. al., Phys. Rev. D41 (1990) 2330 



B. I. Abelev et al. [STAR Collaboration], Ph ys. Rev. C 79, 03 4909 (200977 
STAR collaboration, C. Adlcr et. al., Phys. Rev. Lett. 87 (2001) 112303: PHOBOS collaboration, B. B. Back 
et. al., Phys. Rev. Lett. 85 (2000) 3100 BRA HMS collaboration, I. G. Bearde n et. al, Phys. Lett. B523 (2001) 227 
BRA HMS collaboration, I. G. Bearde n et. al., \Phys. Rev. Lett. 88 (2002) 202301 PHO BDS collaboration; B . B. Back 



et. al., \Phys. Rev. C 65 (2002) 061901: PHENI X collaboration, S . S. Adler et. al., ^Phys. Rev. C71 (2005) 034908| 

R. Baier, A. H. Mueller, D. Schiff and D. T. Son. larXiv:1103.1259l [nucl-th]. 

J. -P. Blaizot, F. Gelis, J. Liao, L. McLerran and R. Venugopalan, arXiv: 1 10 7.5296 1 [hep-ph] . 
B. B. Back et al. [PHOBOS Collaboration], Phys. Rev. C 65, 061901 (2002). 
K. Aamodt et al. [ALICE Collaboration], Phys. Rev. Lett. 106, 032301 (2011). 

B. B. Back, M. D. Baker, D. S. Barton, R. R. Betts, M. Ballintijn, A. A. Bickley, R. Bindel and A. Budzanowski et al, 
Phys. Rev. Lett. 91, 052303 (2003). 

S. Chatrchyan et al [CMS Collaboration], JHEP 1108, 141 (2011). 

B. B. Back et al [PHOBOS Collaboration], Phys. Rev. Lett. 93, 082301 (2004). 
I. Arsene et al [BRAHMS Collaboration], Phys. Rev. Lett. 94, 032301 (2005). 

D. d'Enterria, G. K. .Eyyubova, V. L. Korotkikh, I. P. Lokhtin, S. V. Petrushanko, L. I. Sarycheva and A. M. Snigirev, 
Eur. Phys. J. C 66, 173 (2010). 

W. A. Horowitz and Y. V. Kovchegov, N ucl. Phys. A 849 , 72 (2011). 
G. A. Chirilh, B. -W. Xiao and F. Yuan. larXi^:1112.1061l [hep-ph]. 
A. H. Rezaeian, arXiv:1111.2312 [hep-ph]. 

E. Levin and A. H. Rezaeian, Phys. Rev. D 82, 054003 (2010). 
T. Lappi, Eur. Phys. J. C 71, 1699 (2011). 

J. L. Albacete and A. Dumitru, arXiv:1 011.5161^ [hep-ph]. 

A. Dumitru, D. E. Kharzeev, E. M. Levin and Y. Nara. larXiv:1111.303T] [hep-ph]. 

E. Levin and A. H. Rezaeian, Phys. Rev. D 82, 014022 (2010). 

J. L. Albacete and C. Marquet, Phys. Lett. B 687, 174 (2010). 

J. Jalilian-Marian and A. H. Rezaeian, arXiv;1110.2810 [hep-ph]. 

A. Dumitru, A. Hayashigaki and J. Jalilian-Marian, Nucl. Phys. A 765, 464 (2006). 

T. Altinoluk and A. Kovner, Phys. Rev. D 83, 105004 (2011). 

P. Quiroga-Arias, J. G. Milhano and U. A. Wiedemann, Phys. Rev. C 82, 034903 (2010) 

C. A. Salgado, J. Alvarez-Muniz, F. Arleo, N. Armesto, M. Botje, M. Cacciari, J. Campbell and C. Carli et al, 
afxrv:1105.3919 [hep-ph]. 

G. G. Barnafoldi, J. Barrette, M. Gyulassy, P. Levai and V. Topor Pop. yarXivillll. 36461 [nucl-thj. 



